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Abstract 

We perform a strong-coupling expansion up to the third order of the hopping parameter t for 
the spin-1 Bose-Hubbard model with antiferromagnetic interaction. We examine the validity of the 



expansion by comparing it with the mean- field theory. Moreover, we apply a few extrapolation 



methods up to the infinite order of t to improve our results. 
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I. INTRODUCTION 



Ultracold bosons have been studied extensively since the realization of the Bose-Einstein 
condensation. Here temperature can be approximately zero, which is very difficult to realize 
in conventional experimental systems. In addition to the conventional spinless bosons, spinor 
bosons have been examined as a new bosonic system with the multiple internal degrees of 
freedom The development of the optical lattice systems has further promoted the 

study of ultracold bosons. In particular, the superfluid (SF)-Mott insulator (MI) transition 
has been realized in an optical lattice system |3j-l7fl. 

Theoretically, the optical-lattice system having low boson fillings can generally be de- 
scribed by the Bose-Hubbard (BH) model 0, 0]. In addition, both the MI phases and 



10H28(]. The 



SF-MI transitions of spin-1 bosons have been intensively studied analytically 
ferromagnetically interacting system is essentially similar to spinless bosons, whereas the an- 
tiferromagnetically interacting system has rich physical properties. Several spin phases such 
as the singlet , nematic, and dimeri zed pha ses in the insulating phase have been examined an- 
alytically 
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14j and numerically 
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25] . As for the SF-MI transition, Tsuchiya, Kurihara 



and the present author [116] have used perturbative mean-field approximation (PMFA) 29], 
which expands the free energy in the SF-order parameter, to quantitatively show that the MI 
phase for even boson fillings is considerably stable against the SF phase than that for odd 
boson fillings. This conjecture has been confirmed by the density-matrix-renormalization 
group (DMRG) [22| and quantum-Monte-Carlo simulation (QMC) 24], |25| in one dimension 
(ID). On the other hand, mean-field (MF) studies beyond the perturbation theory 17H19| 
have shown a possible first-order SF-MI transition of the BH model by MF theory for a 
weak antiferromagnetic interaction such as U2/U0 ~ 0.04 corres pon ding to 23 Na. The first- 
order transition in ID has also been shown by a QMC study |24J. On the other hand, if 
the antiferromagnetic interaction is adequately strong, the first-order transition may be ne- 
glected because it occurs when kinetic energy is considerably greater than antiferromagnetic- 
interaction energy near the SF-MI phase boundary 30|. For a second-order transition, the 



strong-coupling expansion method |31] of the kinetic energy, which is based on the Rayleigh- 
Schrodinger perturbation theory [321 ] . is an excellent method for obtaining the phase bound- 
ary. The strong-coupling expansion has been applied in the analysis of the spinless BH 
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-37], extended BH 38], hardcore BH 39], and two-species BH models 40], and the 
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39fl. So far, however, only MF 



results are in excellent agreement with QMC results 
calculations have been performed to analytically study the SF-MI transition of the spin-1 
BH model. 

In this study, we perform a strong-coupling expansion of the spin-1 BH model up to 
the third order of the hopping parameter t. This paper is organized as follows: Section II 
introduces the spin-1 BH model and strong-coupling expansion. Section III gives the results: 
the phase diagrams for two dimensions (2Ds) or three dimensions (3Ds), the critical values of 
t at the peak of the Mott lobes and their dependence on antiferromagnetic interaction, which 
reflect the validity of the expansion, results obtained by several extrapolation techniques up 
to the infinite order of t, and the ID phase diagram. A brief summary of our results and 
discussions are given in Sec. V. 



II. SPIN-1 BOSE-HUBBARD MODEL 

The spin-1 BH model is given by H = H + Hi, 



i,a i,oi,p 

+ \ U z Yl 4*4r F «0 • F lS a iS a i/3 . 

i,a:,/3,7,(5 
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Here, \x and t(> 0) are the chemical potential and the hopping matrix element, respectively. 
Uq (U2) represents the spin-independent (dependent) interaction between bosons. We assume 
Uo and U2 are positive, which corresponds to repulsive and antiferromagnetic interaction. 
The operator (aj Q ) annihilates (creates) a boson at site i with spin-magnetic quantum 
number a = 1, 0, — 1. The number operator at site % is given by rij = Y2a n i& ( n «" = a \ a a ia)- 
Si = p a ia F a paip is a spin operator at site i and F aj 3 represents the spin-1 matrices. 
In this study, we assume a tight binding model with only nearest-neighbor hopping. The 
summation for all the sets of adjacent sites is expressed by For simplicity, we assume 

hypercubic lattices. 



Under the limit of t — > 0, the MI phase exists for arbitrary //; the MI phase also has even no 
bosons per site for Uo(no — 1)— 2U2 < \i < Uqtiq or odd uq boson(s) for UqUq < fi < U$n — 2U2- 
To ensure that the phase diagram has MI phases with odd boson(s) per site, we assume 
Uq > 2U2 > 0. The SF-MI phase boundary can be determined by calculating the energy 
of the MI phase and that of their defect state which has exactly one extra particle or 
hole. Specifically, if E M i{no, /i, t) > (<)E pa,rt (n , fi,t) or E M i{no, > (<)E holc (n , fi,t), 
we determine that the phase is SF (MI), where EMi{no, fi,t) is the energy of the MI state 
and E part (no, /i, t) [E hole (no, /i, t)] is the energy of the defect state with one extra particle 
(hole). The SF-MI phase boundary is determined by 

E Ml (n ,fi,t) = E^ t (n , f i,t) (2) 

or 

E Ml (n ,fi,t) = E holc (n ,fi,t). (3) 



III. STRONG-COUPLING EXPANSION 

By following Ref. jsi|, we employ the Rayleigh-Schrodinger perturbation theory for the 
third order of the hopping parameter t. We start from the unperturbed MI states, define the 
defect state by doping an extra particle or hole to the MI states, and compare the energy of 
the MI state with that of the defect state. 



A. MI states at the zeroth order of the hopping parameter 



The unperturbed MI state with even bosons per site is 



^even = \\ K, 0, 0)fc, 



(4) 



where |rz , 0, 0) ^ implies the boson number iV = Uq, the spin magnitude is S = 0, and the 



spin-magnetic quantum number is S z = at site k. For simplicity, we neg 
MI state including S = 2 states. However, from analytical calculations 
2D variational Monte-Carlo study 
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ect the nematic 
and a recent 



26] , the nematic MI phase for even boson filling may 



exist for a weak \J 2 - 



For the unperturbed MI state with odd bosons per site, we assume a nematic MI state 

^ odd = Y[\n , 1, 0)k- (5) 

k 

Although ^fcn-o = rifc l n o> 1; ±l)/c is degenerate with \I/ dd at t = 0, we can easily find that the 
degeneracy is lifted at finite t, and ^ 00 \d has lower energy at least in the third-order pertur- 
bation of t as expected (Appendix A). This is natural because we assume antiferromagnetic 
interaction. At t — 0, a dimerized state 



\n , 1, 0)i <g> |n , 1, 0)j + |n , 1, <8> |n , 1, 1) 



(6) 



is also degenerate with ^odd, where we assume that no — 1 bosons form a spin singlet in each 
site and extra bosons form dimers between adjacent sites. The dimerized state is considered 
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221-125] . Therefore, the validity of our results 



to be the ground state in ID 
based on ^/ dd for odd boson fillings is basically limited to 2D or 3D systems, although the 
dimerized phase cannot be denied even there. We can easily find that the energy of \&dimer is 
greater than that of \& dd i n the third-order perturbation of t in 2D or 3D (Appendix A). We 
note that \l/ e ven (^odd) is also adopted as the ground-state in PMFA (l6j which is compared 
with our results in the following section. 



B. Defect states 



We define the defect states by doping an extra particle or hole to ^even an d ^odd as 
follows: 



= -j= J2 [/<l"o + 1, 1, 0), ® J] K 0, 0) k 
^ e n = JL J2 [M n o - 1> !> 0)i ® 11 K 0, 0), 

V i kj=-i 



V(>part 

odd 



,7, hole 
W odd 



^^[/,|n o + l,0,0> i ®n|n O) l,0>, 

i kj^i 



(7) 
(8) 
(9) 
(10) 



Here N is the number of lattice sites, and fi is the eigenvector of the hopping matrix 



tij with the highest eigenvalue 



3_l|. In this study, because we assume hypercubic lattices 
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with only nearest-neighbor hopping, fi = 1 and the eigenvalue A = zt, where z = 2d is 
the number of nearest-neighbor sites in the rf-dimensional hypercubic lattice. Although, 
x ^odd < ' h0lC ^ nas 110 other degenerate candidates, h ° lc ' > is degenerate with @P art ( hole ) — 

^7= ^ [|n + (— )1, 1, ±l)j ® rifc^j l n o 5 1) 0)fc] • We find that @P artt ( holc ) the exac t same 
energy with \l>even ( ' hole ' ) up to the third order of t, and we can choose ^cvcn' holc ' ) as the defect 
state. This degeneracy will be lifted when we dope particles or holes with finite density 
because the SF phase will have no magnetization. 



C. Ground-state energies and phase diagrams 

By using ^ eveQ and \l/odd, the energies of the MI state with even or odd boson(s) per site 
are obtained as 



-^MI,cvcn(no) 
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up to the third order of t. On the other hand, the energies of the defect states are obtained 

as 
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up to the third order of t. By equating the right-hand side of Eqs. [TBTfTBl to zero, we obtain 
the SF-MI phase boundary t-fi curve /ieven(^); ^evenW; /^oddW) an< ^ /^oddW) w hich show the 
upper branch (corresponding to particle-doping) or lower branch (corresponding to hole- 
doping) of the phase boundary curve around the Mott phase with even boson fillings, and 
those for odd boson fillings, respectively. 

Figures [IH3] show the phase diagram obtained by the calculation. The MI phase for even 
boson fillings is considerably stable against the SF phase than that for odd boson fillings as 
expected from MF and QMC studies. The area of the MI for even (odd) boson fillings is 
more increased (decreased) for U 2 /Uq = 0.3 than that for U 2 /Uq = 0.15. The critical value 
of t on the phase boundary curve is greater in 2D than that in 3D, because the number of 
nearest-neighbor site z and kinetic energy are smaller at the same value of t. These curves 
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FIG. 1. (Color online) Phase diagram obtained by strong-coupling expansion (Eqs. [T3HT61) for 
U2/U0 = 0.15 in 2D. The solid curves show the results up to the third order of t. Results up to the 
first order (second order) of t are also shown by the blue dashed (green dot-dashed) curve. The red 
dotted curve shows the results obtained by PMFA. 

show the convergence of the strong-coupling expansion from the first to the third order of t, 
and we find convergence is excellent except for that in Fig. [T]for a weak ?7 2 — 0.15£/o in 2D. 

Moreover, we plot the results of PMFA; the plot is exact in infinite dimensions as long as 
the SF-MI transition is of the second order. The area of the Mott lobes obtained by strong- 
coupling expansion are greater than those of PMFA. This difference may reflect the quantum 
fluctuations in the MI phases, which are incorporated (neglected) in strong-coupling expan- 
sion (PMFA). 




D. Consistency in PMFA 



PMFA involves MF decoupling theory using perturbation expansion of the SF order 



parameter. Up to a low-order expansion of the SF order parameter in Ref. [161 ]. the first- 
order SF-MI transition cannot be predicted. However, if the SF-MI transition is of the 
second order, the phase boundary curve is exact in infinite dimensions. Hence, if we expand 
the equation for this phase boundary curve, obtained by the PMFA up to the third order 
of t, it must agree with our proposed strong-coupling expansion in infinite dimensions. The 
results of the expansion of this phase boundary curve obtained by PMFA (Eqs. 30 and 46 
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FIG. 2. (Color online) Same plot as that in Fig. CD for U 2 /U = 0.3 in 2D. 
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FIG. 3. (Color online) Same plot as that in Fig. [T]for U2/U0 = 0.15 in 3D. 
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FIG. 4. (Color online) Same plot as that in Fig. CD for U 2 /U = 0.3 in 3D. 
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for odd MI lobes. These are in exact agreement with our proposed strong-coupling ex- 
pansion. Specifically, we find that the solutions of the equations £?Ml,even(odd) i n o) — 
£ , P art ( hole )(n )def,even(odd) are the same as Eqs. ITTH201 under the limit of z — > oo and t — > 
keeping zt as a constant. 



E. Critical value of t at the peak of the Mott lobe 

In this subsection, we examine the critical value of the hopping parameter t at the peak of 
the Mott lobe (ic)> at which the upper branch of t-fi curve for the phase boundary converges 
with its lower branch. In addition, the dependence of tc on U 2 indicates the range of the 
application of our proposed expansion for the third order of t. 

11 



Figure [5] shows the dependence of tc on U2 for an Mott lobe with an even number of 
bosons (no = 2). The curves of infinite dimensions obtained by PMFA and those obtained 
by the strong-coupling expansion up to the third order of t (see the previous subsection) 
are smoothly increasing functions of U<i- The results obtained by strong-coupling expansion 
up to the third order of t in 2D and 3D show a similar U2 dependence of tc when U2/U0 is 
large. However, they show a uncommon behavior for U2/U0 ~ 0.1, and the curves stop at 
small U2/U0 because we cannot find tc as the upper and lower branches of the t-/i curve 
no longer converge. Such a situation also occurs for a greater n. Because U2 stabilizes the 
MI phase with even boson filling against the SF phase, the results of PMFA and the strong- 
coupling expansion of infinite dimensions are in accordance with physical intuition. A the 
small U2 regime is hazardous for strong-coupling expansion of finite dimensions because a 
few denominators of the expansion have 3U2 which corresponds to the spin-excitation energy 
[E(S = 2) — E(S = 0)] of an intermediate state that appears in perturbative calculation. On 
the other hand, the terms including 3U2 disappear in the denominators of strong-coupling 
expansion in infinite dimensions. 

Figure [6] shows the same plot that in as Fig. [5] for a Mott lobe with an odd number of 
bosons (n = 1). The parameter tc obtained by PMFA is a smoothly decreasing function 
of U2 and becomes zero at U2/U0 = 0.5, at which the Mott lobes with odd boson fillings 
disappear. Moreover, the other curves also show that tc is a decreasing function of U2 
when U2/U0 is large. However, tc is a rapidly increasing function of U2 for U2/U0 < 0.1 
and tc = at U2/U0 = 0. This is also because a few denominators have 3U2 in strong- 
coupling expansion up to the third order of t (here, not only in 2D or 3D but also in infinite 
dimensions). Therefore, tc — > when U2/U0 — » 0, so that the two branches of t-/i curve 
converge. For a greater n , tc cannot be found even for a small U2/U0, same as that for 
the above-mentioned case of even no (figure not shown). These problems can be solved only 
by expanding /i in the infinite orders of t, same as that in PMFA. A similar problem can 
also occur for ~ 0.5 in ID because a few terms in the expansion include Uq — 2[/ 2 i n 

their denominator (see subsection G). The parameter tc may be an increasing (decreasing) 
function of U2 when the Mott lobe has an even (odd) boson filling, and strong-coupling 
expansion can be depended upon only in a large (an intermediate) U2 region. 
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FIG. 5. (Color online) Dependence of tc on U2 for the Mott lobe with an even boson filling uq = 2. 
The blue short-dashed, green dot-dashed, and red long-dashed curves show the results obtained by 
strong coupling expansion up to the third order oft in 2D, 3D, and infinite dimensions, respectively. 
The solid curves show the results obtained by PMFA. 




FIG. 6. (Color online) Same plot as that shown in Fig. [5] for the Mott lobe with an odd boson 
filling no = 1. 

F. Extrapolation methods 

Our expansion up to the third order of t has a few problems. For example, the phase 
boundary curve including the value of tc does not completely converge. From a qualitative 
point of view, expansion does not provide an appropriate scaling form of the phase boundary 
curve near tc- However, in this subsection, we attempt two extrapolation methods toward 



13 



0.15 



0.05 



: u 2 /u =03 








"«0 = 


2 (2D) 


*■* " 




- 
- 




«o = 


2 (3D) 










__ __. .- ~< 




n = 


T(2D) 

— H 










"o = 


1 (3D) 



0.5 1 

1 /order 



FIG. 7. Solid circles show tc obtained by our proposed strong-coupling expansion at the first, 
second, and third orders of t for U2/U0 = 0.3 in 2D and 3D. The data are indicated for no = 1 
Mott lobes and uq = 2 Mott lobes. The horizontal axis shows the inverse of the order of expansion 
(1/order = 1/3, 1/2, 1 for third, second, and first orders, respectively). The dashed lines are least- 
square linear fits for the solid circles. 

an infinite-order theory to improve our phase diagram. 



1. Linear fit of tc 



Critical-point extrapolation, which was proposed in Ref. 31], is a simple method involv- 
ing the use of least-squares fit to the best straight line. 

Figure [7] shows the critical point t c at the peak of the Mott lobe for each order (m) of 
strong-coupling expansion. The data of tc lie approximately on the straight line. The data 
can be extrapolated to the infinite order (1/m — > 0) by using the least-squares fit for the 
best straight line. Specifically, the point crossed by the straight line over a vertical axis is 



the resulting tc by extrapo 
infinite order. The error 



ation to the infinite order. Table [J shows the fitting data to the 
41] of the fitting is very small when Uz/Uq is large, at which strong 
coupling can be more reliable than that in a small U2/U0 regime (which is consistent with 
the dependence of tc on U2 explained in the previous subsection). 
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2. Extrapolation of the phase boundary curves 



Our proposed phase-boundary curve obtained by expansion up to the third order of t has 
a cusp at the peak of the Mott lobe. However, it can be assumed that chemical potential 
has the following power-law scaling behavior near tc, similar to that the spinless BH model 
in 2D or 3D. 

f i = A{t)±B(t)(t c -t) zv . (21) 



31 



38fl. Here, A(t) 



The following fitting method is called as chemical-potential fitting 
a + bt + ct 2 + dt 3 and B(t) ~ a + /3t + ^t 2 + 5t 3 are the regular functions of t. The parameter 
zv is the critical exponent of our model. By using expansion up to the third order of t, we 
immediately determine a,b,c, and d by setting A(t) = |/A part (t) + /i hole (i)]/2. On the other 
hand, we assume that the value of zv to determine B(t). If we assume the scaling behavior 
is the same as that of the spinless BH model or equivalent to that of the (<i+l)-dimensional 
XY model, zv ~ 2/3 for d = 2 and z = 1/2 for d > 2 in the rf-dimensional spin-1 BH model 



42|. By setting 5 = 0, we obtain a, (3, 7, and tc by comparing the Taylor expansion of t 
in B(t)(t c - t) zu with [/i part (t) - /i holo (t)]/2. The values of t c are shown in Table Hand the 
phase-boundary curves obtained by the above fitting are shown in Fig. [8] for 2D and Fig. [9] 
for 3D. 

We can combine these two fitting methods to obtain the phase-boundary curve. Specifi- 
cally, we can use the value of tc obtained by the least-square fit to compare B(t)(tc — t) zu 
with [yU part (t) — /i hole (t)]/2. Here we can include the St 3 term in B(t). The obtained phase- 
boundary curves are also plotted in Figs. [S] and HJ Practically, the phase-boundary curves 
plotted by these two fitting methods are similar, especially in 2D. 

The phase-boundary curves plotted by these two fitting methods are more similar to 
those obtained by PMFA in 3D than those in 2D, as expected. 



G. One dimension 



In ID, the MI phase has very rich spin structures; however, our strong-coupling expansion 
is based on the spin-singlet (spin-nematic) states for even (odd) boson filling. In particular, 
for odd boson fillings, the ground state may be the spin dimerized state in the wide parameter 
regime. Hence, our results cannot be directly applied to ID especially for odd boson fillings. 
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The range of the U2/U0 value, in which we can obtain tc at the peak of the Mott lobe, 
is more limited compared to those in 2D or 3D. For instance, U2/U0 > 0.255 is needed to 
obtain tc for Uq = 2 Mott lobe. 

Nevertheless, we briefly examine the phase diagram because most of the numerical simula- 
tions have been performed for ID so far. Figures ITUl and [TT1 show the phase diagrams obtained 
by our proposed strong-coupling expansion up to the third order of t for U2/U0 = 0.3 and 
U2/U0 = 0.4, respectively. In Fig. [TTJ the upper and lower branches of the no = 1 Mott lobe 
do not converge, and a closed phase-boundary curve cannot be obtained. In addition, Figs 
[TOl and [11] also show the convergence of the strong-coupling expansion from the first to the 
third order of t, and we find that the convergence is not excellent compared with that in 2D 
or 3D (Figs. HJfl). 



For U2/U0 = 0.3, the agreement with the phase diagram plotted by DMRG in Ref. 22 1 



is not excellent but is satisfactory to some extent for the Mott lobe with even boson fillings. 
For example, tc/Uo = 0.327 by our proposed strong-coupling expansion, and tc/Uo ~ 0.47 
by DMRG (the author's reading off Fig. 1 in Ref. [22]). However, for U 2 /U = 0.4, the 
agreement with QMC data, shown in Fig. 1 of Ref. [25] is not satisfactory: tc/Uo = 0.422 
in our proposed strong-coupling expansion and tc/Uo ~ 0.7 in the QMC data. In general, 



TABLE I. List of the critical points t c /Uo. 









Two dimensions 






Three dimensions 




no 


U2/U0 




tc/Ua^ 




(*c/C/ )3rd a 


tc/U b 


tc/Uo" 


1 


0.2 


0.0568 


0.0400 ± 0.0010 


0.0482 


0.0353 


0.0224 ± 0.0020 


0.0269 




0.3 


0.0429 


0.0349 ± 0.0021 


0.0365 


0.0263 


0.0195 ± 0.0001 


0.0202 




0.4 


0.0234 


0.0205 ± 0.0015 


0.0201 


0.0143 


0.0115 ± 0.0004 


0.0110 


2 


0.2 


0.1221 


0.1019 ± 0.0218 


0.1084 


0.0758 


0.0615 ± 0.0078 


0.0607 




0.3 


0.1388 


0.1205 ± 0.0068 


0.1208 


0.0866 


0.0722 ± 0.0021 


0.0683 




0.4 


0.1567 


0.1389 ± 0.0010 


0.1356 


0.0978 


0.0827 ± 0.0009 


0.0766 


3 


0.2 


0.0313 


0.0216 ± 0.0015 


0.0266 


0.0195 


0.0121 ± 0.0016 


0.0149 




0.3 


0.0236 


0.0190 ± 0.0009 


0.0201 


0.0145 


0.0106 ± 0.0001 


0.0111 




0.4 


0.0129 


0.0113 ± 0.0008 


0.0111 


0.0079 


0.0063 ± 0.0002 


0.0061 


4 


0.2 


0.0754 


0.0611 ± 0.0150 


0.0669 


0.0469 


0.0368 ± 0.0057 


0.0375 




0.3 


0.0859 


0.0726 ± 0.0057 


0.0747 


0.0536 


0.0435 ± 0.0021 


0.0422 




0.4 


0.0971 


0.0840 ± 0.0007 


0.0839 


0.0606 


0.0498 ± 0.0003 


0.0474 



a Data obtained by third-order strong-coupling expansion 

b Data obtained by least-squares fit based on the strong-coupling expansion. 

c Data obtained by chemical-potential fit. 



16 



4 



2 





>N 2D 






- "0=4 




** 


- "o= 3 


n =2 











0.05 0.1 



FIG. 8. (Color on line) Phase diagram obtained by strong-coupling expansion up to the third order 
of t (solid curve) for U2/U0 = 0.3 in 2D and its extrapolation to infinite order of t. The blue dashed 
(green dot-dashed) curve shows the chemical-potential fitting without (with) the least-squares fit 
for tc- The blue dashed and green dot-dashed curves are similar to each other. Especially, the two 
curves for the uq = 2 Mott lobe are indistinguishable in the figure. The red dotted curve shows the 
results of PMFA. The area of the Mott lobe is smallest in PMFA. 




FIG. 9. Same plot as that in Fig. [8] for 3D. 



a larger t is required to realize the SF phase for large Uq and/or U2, where the higher-order 
terms of t become more prominent. Therefore, we may have to expand up to fourth- or even 
higher-order terms to reduce the discrepancy. Otherwise, we may have to assume another 
MI phase such as a nematic MI phase, if the MI phase is more stable than the singlet MI 
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FIG. 10. (Color online) ID phase diagram at U2/U0 = 0.3. The solid curves show the results 
obtained by strong-coupling expansion up to the third order of t. Results up to the first order 
(second order) of t are also shown by the blue dashed (green dot-dashed) curve. The red dotted 
curve shows the results obtained by PMFA. 

phase that is assumed in this study. 

Figures [TOl and [TT1 also show the results obtained by PMFA, which may be inaccurate in 
ID. We find that the difference between strong-coupling expansions expansions is very large. 

As mentioned in the previous subsection, we also attempt to extrapolate our results to 
infinite order of t. In the least-squares fit of tc, we cannot improve our results as extensively; 
for instance, t c /U = 0.433±0.030 at U 2 /U = 0.4 for the n = 2 Mott lobe. We also perform 
the chemical-potential fit by assuming the Kosterlitz-Thouless form, as mentioned in Ref. 



3l| . However, we cannot obtain successful fitting results. 



H. Summary and Discussion 

In this study, we have performed strong-coupling expansion of the hopping parameter t to 
analytically obtain the phase diagram. Under the limit of infinite dimensions, the t-\x phase 
boundary curves are consistent with the exact results obtained by PMFA. On the whole, 
the convergence of the phase-boundary curves from the first to the third order is excellent 
in 2D or 3D. The dependence of tc on U2 at the peak of the Mott lobe indicates the validity 
of our expansion. We have attempted to extrapolate our results up to the infinite order of t 
by using the least-squares fitting and chemical-potential fittings of tc- The results are more 
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0.5 
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o 



FIG. 11. (Color online) Similar plot as that of Fig. [TUlfor U2/U0 = 0.4 in ID, except for the dashed 
curves for the no = 1 Mott lobe. For no = 1, the upper and lower branches of the Mott lobe do not 
converge in the third order of t, and we cannot obtain a closed phase-boundary curve. 



equal to those of PMFA in 3D than those in 2D, as expected. In ID, we have compared 
the phase boundary curves obtained by strong-coupling expansion with those obtained by 
numerical simulations. For U2/U0 = 0.3, the agreement with the result obtained by DMRG 
is not excellent but is satisfactory to some extent. 



Our proposed strong-coupling expansion depends on the t = ground state. However, 
the MI phase can be more complicated. We have to consider the dimerized-spin phase for 
odd boson fillings, which can be the ground state in ID. In addition, we should consider the 
nematic spin phase for even boson fillings, which can be the ground state for weak t/ 2 . To 
analytically determine the complete phase diagram, we have to include such spin phases in 
our strong-coupling expansion. This still remains a problem in future. 
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Appendix: Energies of insulating states 



In Sec. Ill A, we have introduced three degenerate ground states at t = when the boson 
filling n Q is odd. 



^odd = Yl l n o, l,0) fc , 

k 

*ferro = K, 1, ±l)fc 

fc 

^dimer = XT (l^O, 1, l)i ® |™0, 1, 



(A.l) 

(A.2) 



-|n , 1, 0}j ® |%, 1, 0)^ + |n , 1, ® |n , 1, 1}/) . 



(A.3) 



By using the Rayleigh-Schrodinger perturbation theory [32], we obtain the energies as fol- 
lows: 



£wn nT 7 (no- l)(n + 4) 2 2n\ + 6n + 7 



y ferro 

N 

-^dimer 



75 



C/n + 4C/ 2 



+ 



15 U Q + U 2 



t 2 fl6(n - l)(n + 4) 2ng + 6n + 7 



+ 



+ 0(t 4 ), 

o K + i)K + 2) 



(A.4) 



27 15 U + iU 2 U + U 2 U - 2U 2 

n -l)(n + 4) 2n ! + 6n + 7 (n + l)(n + 2) 



4 v_u ^ , — + 

£7 + 4[/ 2 U + U 2 U -2U 2 



}+0(t 4 ). (A.5) 



-^Mi.odd i s obtained as specified in Eq. [T21 up to the third order of t. Comparing Eq. IA.4I 
with Eq. [T21 we can easily find that Ef crro > i?Mi,odd for U 2 > at any value of no and z = 2d. 
By comparing Eq. IA. 51 with Eq. [121 we also nn d that -Earner > -^Mi.odd for U 2 > in 2D or 
larger dimensions. On the other hand, in ID, -Edimcr > (<)-EMi,odd f° r U 2 < (>)Uo/6 > 
when no — 1 and -Edimer < E M i,odd f° r U 2 > when n > 3. 



ACKNOWLEDGMENTS 



The author would like to thank Enago (www.enago.jp) for the English language review. 



[1] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H.-J. Miesner, J. Stenger, 

and W. Ketterle, Phys. Rev. Lett. 80, 2027 (1998). 
[2] For reviews, D. M. Stamper-Kurn and M. Ueda, unpublished ([arXiv: 1205. 1888). 
[3] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansen, and I. Bloch, Nature 415, 39 (2002). 

20 



[4] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. Gillen, S. Foelling, L. Pollet, and M. 

Greiner, Science 329, 547 (2010). 
[5] W. S. Bakr, J. I. Gillen, A. Peng, S. Foelling, and M. Greiner, Nature 462, 74 (2009). 
[6] M. M. Endres, M. Clieneau, T. Fukuhara, C. Weitenberg, P. Schaufi, C. Gross, L. Mazza, M.C. 

Banuls, L. Pollet, I. Bloch, and S. Kuhr, Science 334, 200 (2011). 
[7] C. Weitenberg, M. Endres, J. F. Sherson, M. Cheneau, P. Schaufi, T. Fukuhara, I. Bloch, and 

S. Kuhr, Nature 471, 319 (2011). 
[8] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 

(1998). 

[9] M.PA. Fisher, P.B. Weichman, G. Grinstein, and D.S. Fisher, Phys. Rev. B 40, 546 (1989). 
[10] E. Demler and F. Zhou, Phys. Rev. Lett. 88, 163001 (2002). 
[11] S.K. Yip, Phys. Rev. Lett. 90, 250402 (2003). 

[12] A. Imambekov, M. Lukin, and E. Demler, Phys. Rev. A 68, 063602 (2003). 
[13] M. Snoek and F. Zhou, Phys. Rev. B 69, 094410 (2004). 

[14] F. Zhou and M. Snoek, Ann. Phys. 308, 692 (2003): F. Zhou, Europhys. Lett. 63, 505 (2003). 

[15] N. Uesugi and M. Wadati, J. Phys. Soc. Jpn. 72, 1041 (2003). 

[16] S. Tsuchiya, S. Kurihara, and T. Kimura, Phys. Rev. A 70, 043628 (2004). 

[17] T. Kimura, S. Tsuchiya, and S. Kurihara, Phys. Rev. Lett. 94, 110403 (2005). 

[18] K. V. Krutitsky and R. Graham, Phys. Rev. A 70, 063610 (2004). 

[19] K.V. Krutitsky, M. Timmer, and R. Graham, Phys. Rev. A 71, 033623 (2005). 

[20] T. Kimura, S. Tsuchiya, M. Yamashita, and S. Kurihara, J. Phys. Soc. Jpn. 75, 074601 (2006). 

[21] M. Yamashita and M. W. Jack, Phys. Rev. A 76, 023606 (2007). 

[22] M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, and R. Fazio, Phys. Rev. Lett. 95, 240404 
(2005). 

[23] S. Bergkvist, I. P. McCulloch, and A. Rosengren, Phys. Rev. A 74, 053419 (2006). 

[24] G. G. Batrouni, V. G. Rousseau, and R. T. Scalettar, Phys. Rev. Lett. 102, 140402 (2009). 

[25] V. Apaja and O. F. Syljuasen, Phys. Rev. A 74, 035601 (2006). 

[26] Y. Toga, H. Tsuchiura, M. Yamashita, K. Inaba, and H. Yokoyama, J. Phys. Soc. Jpn. 81, 
063001 (2012). 

[27] M. Lacki, S. Paganelli, V. Ahufinger, A. Sanpera, and J. Zakrzewski, Phys. Rev. A 83, 013605 
(2011). 

21 



[28] A. Wagner, A. Nunnenkamp, and C. Bruder, unpublished (larXiv:1207.291ip . 

[29] D. van Oosten, P. van der Straten, and H.T.C. Stoof, Phys. Rev. A 63, 053601 (2001). 

[30] For instance, the first-order SF-MI transition from the Mott lobe with two bosons per site dis- 
appears for U2/U0 > 0.32 in the Gutzwiller MF study as performed in Ref. \v\ (unpublished). 

[31] J. K. Freericks and H. Monien, Phys. Rev. B 53, 2691 (1996); J. K. Freericks and H. Monien, 
Europhys. Lett. 26, 545 (1994). 

[32] See. e.g., C. Cohen- Tannoudji, B. Diu, and F. Laloe, Quantum Mechanics (Wiley-Interscience, 
1992). 

[33] T. D. Kiihner and H. Monien, Phys. Rev. B 58, 14741 (R) (1998). 

[34] P. Buonsante,V. Penna, and A.Vezzani,Phys. Rev. B 70, 184520 (2004). 

[35] K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005). 

[36] J. K. Freericks, H. R. Krishnamurthy, Y. Kato, N. Kawashima, and N. Trivedi, Phys. Rev. A 

79, 053631 (2009). 
[37] V. K. Varma and H. Monien, Phys. Rev. B 84, 195131 (2011). 
[38] M. Iskin and J. K. Freericks, Phys. Rev. A 79, 053634 (2009). 
[39] I. Hen, M. Iskin, and M. Rigol, Phys. Rev. B 81, 064503 (2010). 
[40] M. Iskin, Phys. Rev. A 82, 033630 (2010). 

[41] Here, we adopt standard deviation twice as an error bar because we can use only three data 
to form a straight line, and hence, the accuracy of standard deviation is almost the same as 
standard deviation. [See e.g., Sec. 2.5 of N. C. Barford, Experimental Measurements: Precision, 
Error, and Truth (John Wiley & Sons, 2nd edn., 1985).] 

[42] The critical exponent is determined by phase fluctuation on the fixed-density contour, so it 
can be considered that the only relevant degree of freedom is the SF phase, same as that in 
the spinless BH model. In addition, at U2/U0 —> 0, the critical exponent must be the same 
as that of the spinless BH model because the spin magnetic quantum number of all bosons 
must be the same (there is degeneracy among S z = 0, ±1) to maximize the absolute value of 
kinetic energy, and the spin degree of freedom does not influence the critical behavior. Even if 
U2 becomes finite (but sufficiently small) and degeneracy is lifted, the polar state consisting of 
only S z = bosons (without S z = ±1) remains as the ground state. The critical behavior does 
not usually depend on the value of microscopic parameters, and hence, the critical behavior 
for large U2 is probably the same as XJ% ~ 0. 

22 



